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ABSTRACT 

I 

L  This  report  describes  an  application  of  Augmented  Lagrangian  techniques 
to  the  numerical  solution  of  quasistatic  flow  problems  in  incompressible 
viscoplasticity,  focusing  on  cases  where  the  internal  viscoplastic  dissipation 
potential  is  not  a  differentiable  function  of  the  material  deformation  rate. 
The  stresses  of  elastic  origin  are  neglected,  and  the  variational  formulation 
of  these  problems  is  approximated  via  mixed  finite  elements  of  order  1. 
Convergence  results  are  proved  or  recalled,  both  for  the  finite  element 
approximation  and  for  the  augmented  lagrangian  algorithm.  A  detailed  study  of 
the  local  minimization  problems  which  occur  in  the  augmented  lagrangian 
decomposition  of  the  above  problems  is  also  presented,  together  with  several 
numerical  results.  These  results  were  obtained  using  the  MODULEF  finite 
element  code  on  a  VAX  780  at  the  Mathematics  Research  Center  and  cover 
successively  the  case  of  Norton,  of  Bingham  and  of  Tresca  type  materials. 
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significance  and  explanation 


Augmented  lagrangian  methods,  introduced  around  1970  by  M.  R.  Hestenes 
and  M.  J.  D.  Powell,  are  now  classical  numerical  tools  in  scientific 
computation.  They  take  into  account  the  dual  structure  that  most  problems  in 
continuum  mechanics  do  present,  involving  usually  both  stresses  and 
displacements  (or  velocities),  to  reformulate  them  as  saddle-point  problems, 
which  can  then  be  solved  numerically  by  Uzawa  type  algorithms.  These  methods 
have  already  been  used  in  situations  like  viscoplasticity  by  GLOWINSKI  and 
MAROCCO  [1975]  and  are  described  in  detail  in  FORTIN  and  GLOWINSKI  [1982]. 
Compared  to  previous  publications,  this  report: 

(i)  tries  to  present  a  clean  and  updated  version  of  these  techniques, 

(ii)  uses  a  low  order,  convergent  finite  element  for  the  approximation 
of  incompressible  velocity  fields, 

(iii)  and  studies  in  details  each  local  minimization  problem  which 
appears  during  the  algorithm. 

The  main  mathematical  tool  used  herein  will  be  convex  analysis.  The  goal 
of  this  report  is  to  give  a  comprehensive  presentation  of  all  the  theoretical 
aspects  which  are  behind  the  application  of  augmented  lagrangian  techniques  to 
viscoplasticity  (existence  theory,  approximation,  convergence  of  the 
algorithm,  . . . )  so  that  the  reader  may  be  able  to  implement  these  techniques 
in  any  finite  element  code,  to  obtain  reasonable  numerical  results  with  a 
minimal  experimentation  time,  to  assess  the  validity  of  his  numerical  results 
and  to  judge  the  efficiency  of  his  numerical  technique. 
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The  responsibility  for  the  wording  and  views  expressed  in  the  descriptive 
sumnary  lies  with  MRC,  and  not  with  the  author  of  this  report. 
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NUMERICAL  SOLOTION  OP  VISCOPLASTIC 
PLOW  PROBLEMS  BY  AUGMENTED  LAGRANGIANS 


Patrick  La  Tallac 

1.  INTRODUCTION  AND  FOMULATION  OF  THE  CONTINUOUS  PROBLEMS. 

1.  Introduction.  We  consider  in  this  report  the  problem  of  computing  the  quasiststic 
flows  of  incompressible  viscoplaptic  materials  subjected  to  given  distributions  of 
external  loads .  The  constitutive  law  which  model lies  the  behavior  of  the  considered 
viscoplastic  materials  and  the  configuration  of  the  body  are  supposed  to  be  given.  The 
unknown  is  the  velocity  field  inside  the  body  resulting  from  the  application  of  the 
external  loads. 

The  materials  which  are  involved  in  such  problems  include  freshly  mixed  concrete, 
bitumen,  frozen  soils,  different  types  of  mud,  polymers  at  high  temperature  or  very  hot 
metals-.  These  materials,  whan  subjected  to  external  loads,  fLow  viscously  in  a 
nonreversible  pattern  and  develop  stresses  which  are  mainly  of  viscous  origin.  Most  of 
these  materials  flow  in  an  incompressible  or  nearly  incompressible  way. 

Herein,  to  compute  the  velocity  field  v,  we  use  a  variational  formulation  of  the 
mechanical  problem  (Sec.  1),  which  neglects  the  stresses  of  elastic  origin,  we  discretize 
the  space  of  kinematically  admissible  incompressible  velocity  fields  by  mixed  finite 
elements  of  order  1  (Sec.  2),  and  finally  we  solve  the  resulting  discrete  problem  by 
augmented  lagrangian  techniques  (Sec.  3).  Convergence  results  are  proved  both  for  the 
finite  element  approximation  and  for  the  augmented  lagrangian  algorithm,  and  the  local 
problems  which  appear  in  the  augmented  lagrangian  decomposition  are  studied  in  details  In 
Sec.  4.  Several  numerical  results  are  presented  in  Secs.  5  to  7,  successively  for  Norton, 
Bingham  and  Tresca  type  materials.  The  basic  assumption  in  this  work  is  that  the  internal 
dissipation  potential  associated  to  the  considered  vlscoplastlc  material  is  a  convex. 
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continuous  but  not  necessarily  differentiable  function  of  the  deformation  rate  tensor 
Inside  the  body. 

1.2  The  mechanical  problem.  Depending  whether  we  consider  a  specific  piece  of 
material  with  very  little  motion  or  a  specific  domain  with  incoming  and  outcoming 
materialf  the  configuration  ft  given  in  the  data  of  the  problem  will  correspond  either  to 
the  reference  configuration  or  to  the  present  configuration  of  the  body.  In  this  report, 
we  will  suppose  that  it  corresponds  to  the  reference  configuration  of  the  bodyi  in  other 
words,  we  will  consider  solids  in  small  strains.  The  other  case,  associated  to 
vlscoplaatlc  fluids  flowing  viscously,  is  identical  within  the  replacement  of  the 
lagranglan  coordinates  x  by  the  eulerian  coordinates  X. 

Within  this  convention,  the  unknown  velocity  field  is  determined  by  the  two 
mechanical  equations  below  (PER2YNA  [1966] >: 

constitutive  law  (viscoplastic  incompressible  solid  in  small  strains) 

r  (o(x)  ♦  p  mo  e  a  p^*,  b(»)), 

Tr  (*<▼))  -  0,  *<▼)  -  V2  f  ) , 

virtual  work  theorem  (quasistatic  case) 

f  J  o*(Vw+VwT)/2  dx  «  J  f*w  dx  +  J  g*w  da, 

I  «  0  r2 

l  for  any  w  such  that  m  “  0  on 

These  equations  involve  the  Cauchy  stress  tensor  field  o(x)  and  a  hydrostatic  pressure 
field  p(x).  Here,  the  notations  Vw  and  3  P^  represent  the  gradient  of  the  vector 
field  w  and  the  subgradient  of  the  convex  function  Pj(x,*),  respectively.  In 
addition,  and  denote  the  parts  of  the  boundary  of  (2  where  imposed  velocities 
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«*0  and  imposed  tractions  g  are  applied,  respectively.  Moreover,  for  s  fixed,  the 

internal  dissipation  potential  ( • )  is  a  known  convex  function  of  the  time  derivative 

of  the  linearized  strain  tensor  B.  This  function,  defined  here  over  the  space  of 
NxN 

symmetric  tensors  of  R  with  zero  trace  only  depends  on  the  properties  of  the 
considered  viscoplastic  material  at  point  x.  For  example,  if  we  omit  the  argument  x 
for  simplicity,  Morton  and  Bingham  materials  are  characterized  respectively  by 

(1.1)  Pj  (D)  -  ^  (k/5)P  !  0|P  (Norton) 

(1.2)  Pj(D)  -  u  |d|2  +  /2  g  |d|.  (Bingham). 


1.3  Variational  formulation.  If  we  restrict  the  virtual  work  theorem  to  divergence- 
free  test  functions  w  and  if  we  eliminate  the  Cauchy  stress  tensor  a  using  the 
constitutive  law,  then  the  mechanical  equations  above  correspond,  at  least  formally,  to 
the  variational  problem: 


(1.3) 


Minimize  the  dissipated  energy  rate  J(w)  over  the  set  K  of  incompressible 
kinematically  admissible  velocity  fields. 


where  J  and  k  are  respectively  defined  by 


(1.4) 


(1.5) 


J(i»)  *  ]  V.  (Vj  < ?w+VwT j )  dx  -  j  f*w  dx  -  j  g*»  da, 
n  n  r2 


K  “  {w  e  w',P(fl),  div  w  •  0,  w  •  u  on  T.} 

O  1 


This  variational  problem  is  well-posed  and  we  have: 

EXISTENCE  THEOREM :  Let  S5  be  open  bounded  connected  in  R  (N-2  or  3)  with 
Lipschltz  continuous  boundary  T .  We  suppose  that  the  interior  of  T  is  not  empty, 
that  la  the  trace  of  a  function  of  n''p(f!),  and  satisfies 


« 
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I  a  •  v  d«  •  0 

r  ° 

whenever  r  -  r.  We  assume  moreover  that  the  external  body  forces  f  and  surface 

D*  p*  *  * 

tractions  g  are  respectively  In  L  (8)  and  I.  (T  )  (pp  -  P  +  p  ),  and  that  the 
convex  Internal  dissipation  potential  P1  satisfies! 


(1.6)  c1  |d|P  <  P^D)  <  c2  +  c3  |d|p. 


almost  everywhere  In  8  for  any  symmetric,  N  x  N  matrix  D  with  sero  trace. 

1  <  p  <  +-. 

Then,  there  exists  a  velocity  field  v  which  minimizes  the  dissipated  energy  rate 
J(w)  over  the  set  K  of  kinematically  admissible  velocity  fields.  This  solution  Is 
unique  If  P1  Is  strictly  convex.  Moreover,  for  each  mlnlmlzer  v,  there  exists  a 
devlatorlc  stress  tensor  field  c^  in  (Lp  (8))Nx,,,  a  hydrostatic  pressure  field 

p* 

P  in  If  (8)  which  satisfy  the  weak  equilibrium  equations  and  constitutive  laws; 


I  (o  -  pld)'D(v)  dx  *  ]  f*w  dx  +  J  da  ,  w  w  e  V, 

0  f2 


(1.7)  / 


0D  e  (i(v))  a.e. 


in  8, 


|^V  -  {w  e  W1,p(8),  w  -  0  on  ri  >*. 


Proof ;  The  proof  of  this  result  Is  very  classical  in  convex  analysis.  The  existence 
of  a  solution  ▼  Involves  the  Weierstrass  theorem,  its  characterisation  by  (1.7)  uses 
duality  arguments  and  the  closed  range  theorem. 

First,  from  (1.6)  and  from  the  Korn's  inequality  on  V  (GEYMONAT,  SUQUET) 1983) ) , 

J  satisfies 

C,lwlP  -  (Ifl  ,+lgl  ) Iwl  <  J(w)  <  C  (mes(8) )  +  C,lwl  P 
1  1,p  p*  p*  1,p  2  3  1  ,p 

+  (Ifl  ,  +  Igl  ,)M. 

P*  P*  l.p 
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for  any  w  in  K.  Therefore,  J  is  coercive,  convex  (strictly  convex  if  V^(')  is)  and 
continuous  on  K  for  the  w',p(£5)  topology.  It  is  thus  weakly  lower  semicontinuous  on 
K.  In  addition,  K,  defined  as  the  Kernel  of  the  linear  application 
w  ♦  {dlv  w,  w  -  uq  |  p  } ,  is  convex  and  closed  in  -  Since  K  is  also  not  empty, 

applying  the  Weierstrass  theorem,  there  exists  a  minimiter  v  of  J(»)  over  K,  which  is 
unique  if  J  is  strictly  convex. 

To  further  characterize  such  a  minimiter  v,  we  now  introduce 


»  {  »  e  w  «  0  on  ri#  div  v  “  0}, 


Y  -  {D  e  (LP(n))NxN,  DT  «  D,  Tr  (O)  -  0  a.e.  in  15}, 


9(w,D)  »  l  V .(D(v+sr)  -  D)  d x  -  j  f *  ( v+w)  dx  -  J  g  •  (v+w)  da  . 

a  1  a  r 

2 

Above  D(v+»)  represents  the  tensor  D(rt«r)  -VjfVIv+w)  *  V(s*h»)T)  and  Y  is  in  duality 


with  the  space 


Y*  -  {t  e  (Lp"((5))NxN,  tT  -  T,  Tr  (T)  »  0}, 


through  the  duality  pairing 

N 

<  T,  D  >  -  i  T'D  Ox  -  J  l  T.  D  dx  . 

n  n  i,j-i  13  13 

Obviously,  from  (1.6),  *(',•)  takes  on  finite  values  and  is  real,  convex  and  continuous 

on  XxY.  Moreover,  since  ▼  minimizes  J  over  K,  0  is  a  solution  of  the  primal 
problem:  Minimize  9(w,0)  on  X,  From  a  basic  theorem  of  convex  analysis  (EKELAND- 

TIMAM  (19?6,  p  52-53]),  this  implies  that  the  dual  problem:  Maximize  -  9*(0,t)  over 


Y*  has  a  solution  (-o^)  which  satisfies 
(1.8)  {o,  -  d0)  e  39(0,0) 

that  is  <  -  0D,  D  >  <  *(w,D)  -  9(0,0),  V  (*,D)  e  XxY. 

Writing  (1.8)  successively  for  (w,D)  *  Cw,0(v)}  and  {w,D}  •  {0,-B},  we  obtain 


(1.9)  L(v)  -  ]  a  *D(w)  dx  -  1  (•«  dx  -  J  g*w  da  -  0,  V  v  e  X, 

n  D  n  r2 

(1.10)  j  a  *H  dx  <  j  {p  (O(v)  +  H)  -  V.  (D(v)  ) )  dx,  V  B  e  Y. 

fl  0 


But  (1.10)  can  only  hold  If  (EKELAND-TEMAM  [1976  p  21,  p  271)) 

(1.11)  oD  e  aPt(D(v)),  a.e.  in  fl. 

Now,  to  obtain  (1.7)  out  of  (1.9),  (1.11),  it  is  sufficient  to  observe  that  the 
divergence  operator  is  a  continuous  surjection  from  V  onto 
Lp(fl)  (or  onto  Lp(fl)/R  if  I*  »  T).  Therefore,  from  the  closed  range  theorem,  its 
transpose  is  a  continuous  homemorphism  from  L  (0)  onto  the  orthogonal  of  its  Kernel 
in  V*,  that  is  onto  X*.  Since,  from  (1.9)  L( • )  is  an  element  of  X*,  there  exists 

P* 

then  an  element  p  in  L  (fl)  such  that 
L(w)  ■»  <  p,  div  w  >,  V  w  in  V, 
and  our  proof  is  complete. 


REMARK  1 . 1 i  We  are  not  supposing  here  any  differentiability  of  the  Internal 
dissipation  potential  D ^  < • ) .  The  numerical  techniques  to  be  used  later  will  have  to  be 
able  to  handle  such  a  lack  of  differentiability.  q 

REMARK  1.2;  Even  though  the  argument  x  has  been  omitted  in  the  potential  P1  for 
simplicity,  the  whole  theory  developed  in  this  report  applies  for  potentials  which  are 
measurable  functions  of  x  on  (2.  q 

RBI  ARK  1.3:  For  Norton  and  for  Bingham  materials,  the  internal  dissipation  potential 
is  strictly  convex  and  satisfies  (1.7)  with 


C,  -  0,  C  «  C,  -  -  (k/2)p  ,  (Norton) 
2  1  3  p 


c,  -  I), 


:2  *  g,  Cj  »  (/2  9  +  p  )  , 
-6- 


p  “  2  (Bingham) 


But  (1.7)  is  still  valid,  and  therefore  the  above  existence  theorem  still  applies  for 


materials  associated  to  non-strictly  convex  and  non-dif ferentiable  potentials  such  as 

0  <D)  -  -(k/2)p  sup  (|d  -d  |)p, 

P  i,j  3 

where  D^  are  the  eigenvalues  of  the  deformation  rate  tensor  D.  This  corresponds  to 
Tresca's  type  viscoplasticity. 

2.  THE  DISCRETE  PROBLEMS. 

2.1  The  discrete  spaces.  The  approximation  of  the  set  K  of  kinematically 
admissible  velocity  fields,  which  is  needed  for  the  numerical  solution  of  the  variational 
problem  (1.4),  can  not  be  achieved  by  the  basic  finite  element  spaces  used  in  general. 

For  example,  the  space  of  divergence-free  functions  whose  restriction  to  each  triangle 
(tetrahedron  if  N  »  3)  of  a  given  regular  triangulation  of  0  is  a  first  degree 
polynomial  may  only  approximate  a  small  part  of  the  space  of  divergence-free  elements  of 
w1,p(n).  Therefore,  it  is  a  very  inappropriate  finite  dimensional  approximation  of  the 
set  K  of  kinematically  admissible  incompressible  velocity  fields.  To  obtain  a 
satisfactory  approximation  of  K,  the  set  of  approximate  test  functions  must  be  enriched 
and  the  incompressibility  constraint  must  be  weakened. 

As  pointed  out  in  BREZZI  [1974]  and  summarized  in  GIRAULT-RAVIART  [1979]  in  their 
study  of  the  Stokes  problem,  a  good  approximation  of  K  is  obtained  as  follows: 

(i)  we  first  decompose  the  domain  S)  into  a  regular  triangulation  of  Nh 

polygons  (N«2)  or  polyhedrons  (N»3)  which  satisfy  the  classical  assembly  conditions 
described  in  CIARLET  [1978  p  51 ] ; 

(ii)  we  then  define  the  space  Vh  of  approximate  test  functions  by 


wherr  )  is  a  given  finite  dimensional  space  of  continuous  interpolating  functions 

defined  on 

(iii)  in  addition,  we  introduce  an  appropriate  finite  element  space  P^,  included 
00 

in  L  (SI)  and  which  satisfies  the  so-called  BREZZI  (or  inf-sup)  'ondition: 


where  8  is  independent  of  the  diameter  h  of  the  triangulation  T^,  and  where  p  is 
the  exponent  which  appears  in  the  definition  of  K(pp*»p+p*) > 

(iv)  we  finally  approximate  K  by: 


(2.3)  1^  -  ,  (»h“io)  e  Vh,  J  div  dx  -  0,  V  e  P^} 

Briefly  speaking,  this  construction  of  amounts  to  impose  the  incompressibility 

constraint  in  an  averaging  sense  only.  In  that  way,  more  elements  of  Vh  can  satisfy 
this  constraint  and  the  set  is  bigger.  It  can  then  better  approximate  K. 

The  choice  of  the  polyhedrons  of  the  interpolating  space  )  and  of  the 

space  Ph  of  approximate  pressures  is  free,  provided  that  the  BREZZI  condition  (2.2)  is 
satisfied.  In  this  report,  we  will  use  triangles  (respectively  tetrahedrons  if  R  *  3)  as 
polygons  and  define  P^ISl^)  and  P^  by 


(2.4) 

(2.5) 


W  -  {we  c°<nJt) , 
Ph  -  {q  e  C°(C),  q|^ 


w  4  e  r An.), 
1  4 


e  p^Oj),  v  i 


Vi-  1,2N}, 

V' 
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4 


where  is  space  of  first  order  polynomials  defined  over  and  where  (fl^ ) 

U 

are  the  2  triangles 

x  position  of  the  degrees  of 

freedom  for  the  pressures 

o  position  of  the  degrees  of 

freedom  for  velocities 

Figure  2.  1  Decomposition  of  a  triangle 
in  four  equal  subtriangles 


(pressures) 


2 

Figure  2.3  Triangulation 
(velocities) 
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(respectively  tetrahedrons)  Included  in  fig  which  are  obtained  by  joining  together  the 


midsides  of  every  edge  of  fig  With  that  definition  of  t*>*  BPace  vh  of 

approximate  test  functions  is  now  simply  the  apace  of  continuous  vector  functions  with 
zero  trace  on  I"  and  whose  restriction  to  each  triangle  (respectively  tetrahedron)  of 
is  a  first  degree  polynomial ,  Tjj  being  the  triangulation  obtained  by  dividing  each 
triangle  (respectively  tetrehedron)  of  Th  into  four  equal  subtriangles  (respectively 
eight  subtetrahedrons ) •  As  for  the  space  ^  of  approximate  pressures ,  it  becomes  the 
space  of  continuous  scalar  functions  whose  restriction  to  each  triangle  (respeclvely 
tetrahedron)  of  7^  is  a  first  degree  polynomial  (see  GLOWINSKI  (1984)  for  more  details 
on  those  discrete  spaces) • 

The  above  choice  of  approximate  spaces  ((2.1),  (2.3),  (2.4),  (2.5))  is  far  from  being 

the  only  possible  one  but  it  satisfies  the  BREZZI  condition  (2.2)  and  leads  to  a  very 

convenient  approximate  augmented  lagrangian  decomposition  of  our  viscous  flow  problem 
(1.3).  Moreover,  it  uses  low  order  finite  elements,  which  is  adviseable  in  nonlinear 

problems  where  little  regularity  is  to  be  expected.  Finally,  the  sets  (K^)  constructed  by 

(2.1),  (2.3),  (2.4)  and  (2.5)  form  a  converging  sequence  of  finite  dimensional 
approximation  of  K  and  we  have  ( BERCOVIER-PIRONNEAU  [1977]): 


(2.6) 


V  w  e  K,  lim  1  Inf  Iw 
h*°  -h^h 


Vl.p  f 


0  , 


RP4ARK  2. 1 :  When  the  maximal  diameter  h  of  the  triangulation  goes  to  zero,  we  also 
have  (CIARLET  [1978]): 


V  C  V  O»’,q(0),  V  1  <  q  <  +  dim  V  <  +-» 


(2.7) 


Vw  e  V,  lira  {inf  Iw  -  w.  I  } 

w  a\j  h  1,p 


0: 


h*0  w^ev 
n  n 
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(2.8) 


P  C  LH(fl),  V  1  <  q  <  ♦  dim  P  <  +-; 


¥q  e  LP*(fl),  lim  | Inf  Iq  -  a  I  l  -  0. 

h*°  V*h  °'P 


But,  since  we  are  imposing  the  incoapressibility  constraint  in  an  averaging  sense  only, 

*h  is  not  included  in  1C.  q 

2.2  The  discrete  problems .  The  approximate  incompressible  viscous  flow  problem  is 
simply  obtained  by  replaing  K  by  in  (1.3).  But,  since  is  not  included  in  K, 

we  first  have  to  extend  the  internal  dissipation  potential  D^.),  initially  defined  as  a 
convex  continuous  coercive  function  on  the  space  of  symmetric  N  x  N  real  matrices  with 
zero  trace,  to  a  convex  continuous  coercive  function  D®( . )  defined  on  the  whole  space  of 
symmetric  N  x  M  real  matrices.  This  extension  must  be  convex  and  satisfy 

ID®  (D)  -  P^D).  v  D  €  with  Tr  D  -  0, 

D®(DHiM>  >  D^d),  «q  «  I,  V  D  e  l£XN  with  Tr  D-0| 

C,|o|P  <  o*  (D>  <  C2  +  C3  |d|P,  V  D  e  »^N,  a.e.  in  Q  . 

In  other  words,  the  extension  D® ( •  )  of  < • )  coincides  with  P^(*)  on  the  space 

of  symmetric  matrices  with  zero  trace,  penalizes  the  slightly  compressible  velocity  fields 
NxN 

and  extends  to  the  coercivity  and  the  continuity  of  the  function  P  ^ ( •  ) .  The 

introduction  of  D*(«)  does  not  affect  the  solutions  of  the  viscous  flow  problems  (1.3) 
since  D®(»)  and  P^(')  coincide  for  symawtric  matrices  with  zero  trace;  it  only 
provides  a  mathematical  tool  for  calculating  the  dissipation  potential  for  compressible 
velocity  fields  and  therefore  enables  us  to  compute  reasonable  nearly  incompressible 
finite  dimensional  approximations  of  the  solutions  v  of  (1.3).  The  introduction  of  such 
extensions  P®(*)  which  satisfy  (2.9)  is  easy.  For  example,  for  Norton  and  Bingham 
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material*,  the  expressions  of  given  in  (1,1)  and  (1.2)  define  auch  extensions.  In 

other  caaea,  one  can  take 

t?“(D)  -  (^(D  -  Tr(O)Zd)  c,  |tt  (D)|P, 

being  careful  not  to  chooae  too  big  value*  for  Cj,  in  order  to  avoid  'locking*  phenomena. 


Once  thi*  extenalon  defined,  the  dlacrete  variational  formulation  of  our 
lncompraaalble  vlacoua  flow  problema  (1.3)  ia 


(2.10) 


Minimise  the  dlaalpated  energy  rate  J(w^)  over  the  aet  of  approximate 

kinematically  admissible  velocity  field*. 


where  ia  the  aubaet  which  i*  defined  by  (2.1),  (2.3),  (2.4),  (2.5)  and  where  the 

dissipated  energy  rate  J(  • ) ,  given  by  (1.4),  is  extended  to  a  function  from  into 

H  by  replacing  the  internal  dissipation  potential  ( • )  by  ita  extenalon  P*( • ) 
introduced  in  (2.9). 

THKOMM  2. 1 1  Under  the  assumptions  of  Theorem  1.1,  for  any  fixed  h,  the  discrete 
Incompressible  viscous  flow  problem  (2.10)  has  a  solution  vh.  Moreover,  any  solution 
wh  o £  (2.10)  is  associated  to  an  approximate  stress  tensor  field  and  to  an 

approximate  pressure  field  p^  in  Ph,  such  that 


f  1<Jh+»h  M)  e  3  P®(«^)  on  !5,  ^  -  V2<*V*r*). 


(2.11) 


V<W/2  *  "  ■'/'■h  +  >a  r\  da,  V  ^  e  .h 


In  (2.11),  the  approximate  devlatorlc  atreaa  tensor  (oh+p^  M)  belong*  to  the  apace 
*h  of  piecewise  constant  symmetric  matrix  fields  i 


(2.12) 


{  n  id  + 


_NxN  T 
R  ,Dh 


°h'  *V,  , 


e  <p  (ah)NjcN, 

O  l 


Vi 


1.2",  V  t 


1. 


V- 


Proof :  this  theorem  is  the  discrete  equivalent  of  the  existence  Theorem  1.1.  Dp  to 
(1.9),  (1.10)  its  proof  is  identical,  after  replacement  of  by  {’'(•I,  of  K  by 

K^,  of  Y  by  Yh  and  of  X  by 

(2.U,  xh-{«hevh,  /fl%div^-o.  v^ePh>. 

Now,  since  Yh  is  made  of  piecewise  constant  matrix  fields,  (1.10)  will  also  give 

-  (o  >  6  *  P,  (D(vh>)  a.e.  in  0. 

To  finish  the  proof  of  (2.11),  we  introduce  the  operator  B  frost  into  PJ  defined 

by 

(2.14)  <  bv  gh  >  -  Jo  ^  div  ^  dx.  V  ^  e  Ph,  %  e  V 

whose  Kernel  is  X}],  by  definition.  Frost  the  BREZZI  inequality  (2.2),  (see  for  example 
GIRAULT-FAVIART ( 1979,  p  41]),  this  operator  is  a  continuous  surjection  from  onto 

P{J.  Using  the  closed  range  theorem,  its  transpose  is  a  one-to-one  hosMaxorphiam  frost 

*  * 

Ph  onto  the  orthogonal  of  ^  in  Vh>  But,  frost  (1.9),  the  element  Lh<*>  of  vh 

defined  by 

(2.15)  VV  -  J  (oD)h  •  D(^)dx  -  da  -  J  fM^  dx 

belongs  to  this  orthogonal  subspace.  Therefore,  there  exists  a  unique  pressure  field 
p^  in  Ph  such  that 

l^i  (e^|)  »  <BTp^,^  >  »  j  div  dx,  1^  e  V. 
which  is  exactly  (2.11).  Q 

2.3  Convergence  result.  In  order  to  check  that  the  discrete  problem  (2.10)  is  a 
good  approximation  of  the  continuous  viscous  flow  problem  (1.3)  when  the  maximal 
diameter  h  of  the  triangulation  Th  goes  to  tero,  one  must  study  the  behavior  of  the 
sequence  (v^)  of  solutions  of  (2.10)  when  h  goes  to  zero.  We  will  prove  in  this 
paragraph  that  ( v^)  converges  weakly  towards  solutions  ▼  of  the  continuous  problem  and 
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1 


I 


that  the  dissipated  energy  rate  J(v  )  converge#  toward#  J(v)»  Moreover,  under 

h 

additional  uniform  convexity  assumptions,  such  as  those  satisfied  by  Norton  or  by  Bingham 

materials,  there  is  strong  convergence  of  (v  )  towards  v  in  B1  'P(fl ) .  The  next 

n 

theorem  summarises  these  convergence  properties,  denoting  by  q  the  maximum  of  p  (p  is 

the  exponent  introduced  in  (1.6))  and  2  and  by  Yp  the  space 

Y  -  {D  e  (LP(fl))NxM,  DT  -  D}. 

P 


THEOREM  2.2  a  Under  the  assumptions  of  Theorem  1.1,  the  sequence  ( vh )  of  solutions 
of  the  discrete  problem  (2.10)  decomposes  Itself  into  subsequences ,  each  of 
them  converging  weakly  in  **'P(fl)  towards  a  solution  w  of  the  continuous 
incompressible  viscous  flow  problem  (1.3),  when  h  goes  to  zero.  The  dissipated  energy 
rate  J(v  )  also  converges  towards  J(v).  Moreover,  if  the  extended  Internal 
dissipation  potential  p®  is  of  the  form 
(2.16)  p*(D)  -  V,  (»>  “  G  (D)  +  G,(D), 

with  Gi  convex  and  bounded  below,  g0  convex,  differentiable  and  satisfying 


(ID!  ♦  IGI  )q_p  I 
o,p  o,p 


‘35*  <c) 


!fo 

3D 


( D) )  *  ( G-D) dx  >  C  IO-Olq  , 
o  o.p 


V{D,C)  e  y2, 
p 


(2.17)/ 


3G  3G 

J(3T<c> 


3D 


t,P 

then  the  whole  sequence  (v^)  converges  strongly  in  M  (0)  towards  the  unique 
solution  v  of  the  continuous  problem  (1.3). 

Proof i  The  proof  is  an  Immediate  generalization  of  the  techniques  used  by 
GLOWINSKI-LIONS-TREMOLIERES  [1981,  p  361]  in  their  study  of  Bingham  fluids.  It  requires 
three  steps,  step  1  (Vj,)  is  bounded  uniformly  in  h.  Let  v  be  a  solution  of  the 
continuous  problem  (1.3)  and  let  Sj,  be  the  element  of  such  that 
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* 

i  < 


(2.18) 


,v  -  Vi.p 


*»  Inf  It 


Vl.p 


1  -D 

Proa  (2.6)  (a^)  strongly  converges  towards  v  in  w  (ft)  as  h  goes  to  zero.  Since, 

1(D 

by  extension,  J(*)  is  continuous  on  H  (ft),  this  implies  that,  for  h  sufficiently 
small,  we  have 

(2.19)  J(s^)  <  J(v)  +  1  . 

But  since  is  a  solution  of  (2.10),  we  get 

J(vh)  <  J(s^)  <  J(r)  *  1  . 

Prcei  the  convexity  of  J(*),  this  implies 

(2.20)  J<<*  -*  )/2)  <  (J(-i  )  +  J(v)  +  1)/2  -  C,, 

no  o  b 

where,  as  usual,  the  notation  represents  strictly  positive  numbers  independent  of  x 

and  h.  From  (2.9),  P*  is  coercive,  thus  (2.20)  implies 


l<7(W 


+  7  ( 


VVT)/21 


<  c. 


(If l  +  lgl )lv  -u 
h  o 


1.P 


which,  from  the  Korn's  inequality  and  since  p  is  strictly  greater  than  1,  can  only  hold 
if 


Step  2  weak  convergence  of  (vh).  Since  the  sequence  (»h)  is  uniformly  bounded  in 
,P(ft ) ,  it  decomposes  itself  into  subsequences,  each  of  them  weakly  converging  in 
W1,p(ft).  We  still  denote  by  (vh)  such  a  subsequence  and  denote  by  ▼  its  weak 
limit.  Moreover,  let  v  be  a  solution  of  (1.3)  and  let  (Sj,)  be  the  sequence  of 
elements  of  K^  defined  by  (2.18).  Since  vh  minimizes  J  over  1^,  we  have 

(2.21)  J(vh)  <  Jt^)  ,  W  h  . 

Going  to  the  limit  in  (2.21)  as  h  goes  to  zero,  and  using  the  weak  lower  semicontinuity 
of  J  on  the  left-hand  side,  the  strong  continuity  of  J  on  the  right-hand  side,  we 
obtain 
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(2.22) 


J(v)  <  lim  inf  J(v.)  <  lim  sup  J(v)  <  lim  J(m  ) 
n  h  h 

h*0 


J(v). 


On  the  other  hand,  let  q  be  any  element  of  1/  (0)  and  let  (^)  be  the  sequence 
of  elements  of  Pj,  which  strongly  approximates  q  in  L  (0) .  Since  *h  belongs  to 
K^,,  we  have 

(2.23)  J  a  div  v  dx  »  I  a  div(  v-v.  )dx  ,  V  h. 
a  n  ^  h 

Going  to  the  limit  in  (2.23)  as  h  goes  to  zero  and  using  the  strong  convergence  of 
(q^,)  and  the  weak  convergence  of  ( v^ )  yields 

j  q  div  v  dz  ■  0,  Vq  e  LP  (d). 

n 

Moreover,  from  the  weak  continuity  of  the  trace  operator,  v  -  u  has  zero  trace  on 

o 

So,  finally,  ▼  belongs  to  K.  But  v  minimizes  J  over  K,  therefore 

J(v)  <  J(v). 

Combined  with  (2.22),  this  implies  that  J(T)  is  equal  to  J(v)  and  that  all 
inequalities  in  (2.22)  are  equalities,  therefore  ▼  is  also  a  solution  of  (1.3)  and  the 
whole  sequence  J(vh)  converges  towards  J(v). 


Step  3  Strong  convergence  of  ( v^) .  From  now  on,  we  suppose  that  the  extended 
internal  dissipation  potential  p®  satisfies  (2.16)  and  (2.17).  From  (2.17),  G0  is 
strictly  convex  on  tf6*,  therefore,  by  addition,  p®  is  strictly  convex.  So  is  its 
restriction  P1  on  the  space  of  symmetric  matrices  with  zero  trace.  From  Theorem  1.1, 
the  solution  v  of  (1.3)  is  then  unique.  Thus  the  only  possible  weak  cluster  point  for 
the  sequence  (v^)  is  v  and  the  whole  sequence  (Vj,)  of  solutions  of  (2.10)  converges 
weakly  towards  v  in  W1,p(fl). 

To  prove  its  strong  convergence,  we  first  write  the  discrete  weak  equilibrium 
equations  (2.11)  and  the  continuous  weak  equilibrivsn  equations  (1.7)  for  the  test 
function  ”  *h  "  Th  where  le  the  element  of  1^  defined  in  (2.18).  This 

gives 


I 


m 
1 


A 
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(2*24) 


f  jA  v  D(  w*  ■  L  **(  w*  + 1  «*(w  da- 


1  o  •  D(*h-rh)dx  -  J  f  •  (*h-^h)fc  +  J  9  •  <W  0*. 


(o  +  p  Id)  e  3  £>,(*),  *  -  V2  (Vv+VvT), 


V. 


(0h  +^M)  e9  P*<V'  \  -V2(VV’<). 


where  the  notation  D(w)  represents  as  usual  the  symmetric  component 
T  NxN 

(Vtr+Vw  )/2  of  the  matrix  Vw  in  R  .  By  definition  of  the  subgradient,  the  third 
line  of  (2.24)  is  equivalent  to 

(2.25)  (a+p  M)  •  H  <  P  (i+«)  -  V  (■) ,  V  ■  e  RNxN  with  Tr(«)  -  0. 

II  8 

Since  (a+p  Id)  has  zero  trace  and  since  p®  is  an  extension  of  P,  which  satisfies 
(2.9),  (2.24)  yields 

(a+p  U1MB+4  M)  -  (a+p  M)  •  ■  <  P  (■+■)  -  P^i) 

<  P®  (i+P+q  Id)  -  P®  (*),  v  q  e  *, 


or,  in  other  words 
(2,26)  (a+p  U)  e  9  p*  (K)  . 
Now  setting 


9G 


9G 


t  -  o  +  p  Id  -  —  (i),  th  -  ah  +  ph  Id  -  (*,)  . 


we  obtain  by  substractlon  from  (2.24)  and  (2.26) 
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(2.27) 


3  G  36 

lt-V  •  D(W  *  ■  'Q  ‘if  <V  *  3D2  (in  *  *W 


+  L  H1  div  ‘v-h 


T  )  da  , 


t  e  3q  ,<*>  .  \  e 


But,  since  t  and  t^  are  subdifferentials  of  G,,  we  have  by  definition 


(2.28) 


J  t  •  W»h-r)  da  <  J  {Gf(^,)-6f(*)> 


Ja  %  •  °<w  *  * 


A  suitable  combination  of  (2.27)  and  (2.28)  yields 


36  3G  36  36 

(2.29,  (0(^)1  -  ^  (i^-ma^d,  <  »(«,)»  -  ^(in-D  <vV  <*, 


♦  J  (p-s^)  div  dx  +  J  +  G^w^l)  ■ 


Both  *h  and  a^  are  elements  of  Xj,,  so  we  can  replace  in  (2.29) 


J  (p-Ph)  div  da  bY  J  (P*^ >  div  da  , 

where  qh  is  the  element  of  Ph  which  approximates  p  in  L  ( ft ) .  Once  this 

replacement  done,  we  have  from  (2.29),  (2.17)  and  the  Korn's  inequality 


,2-30)  «,."V,/V,.p»”  'WlP  «  c„,Vh,,,P'^U+1,^,i,P  +  ,v,i.p,q'2 


+  ’P^h'cpv^-^h'l.p  +  C12|  i  <f»(’-V  +  G,(D(^),  -  G,(«)}  da|  . 
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Since  by  construction  and  q^  converge  strongly  respectively  towards  v  and  p  in 

W1,P<n)  and  LP  (fl),  since  from  Step  2  <vh)  is  uniformly  bounded  in  W1,P(f!)  and 
since,  from  (2.9),  the  integral  of  G1  is  continuous  on  Y  ,  the  right-hand  side  of 
(2.30)  converge  towards  O  when  h  goes  to  zero.  Therefore  1*^  -  v^l ^  must  also 

converge  to  zero,  and  from  the  triangular  inequality,  the  sequence  strongly  converges 

1  #D 

towards  v  in  W  (ft)  when  h  goes  to  zero.  _ 


REMARK  2.2:  Ttoree  facts  are  crucial  in  our  proof  of  convergence:  the  existence  of  an 
approximate  pressure  the  existence  of  a  sequence  (*^)  of  elements  of 

approximating  v  and  the  existence  of  a  sequence  (q^)  in  approximating  any 

p* 

element  q  of  1.  (SI).  Although  not  necessary  the  BREZZI  condition  (2.2)  is  a  basic 

tool  for  proving  the  first  two  facts.  I-. 


REMARK  2.3:  Norton  and  Bingham  materials  satisfy  the  uniform  convexity  assumptions 

(2.16)  and  (2.17)  (SCHEURETH  {1977J,  GLOVZNSK1-HAROCCO  f  7975) )  and  therefore,  strong 

convergence  can  be  proved  in  both  cases.  Moreover,  the  speed  of  convergence  of  ( vh ) 

towards  v  can  easily  be  estimated  by  (2.30)  as  a  function  of  the  quantity 

[Vq 

using  the  Idpschitz  continuity  of  G1  and  the  Identity 


,9  Kqw 

ab  <  —  +  — r 

q  q* 


Inf  1w  —vl  I 
_  h  1,p) 

K®h  / 


In  addition,  since  the  dissipation  potential  is  continuously  differentiable  for 
Norton  materials,  the  strong  convergence  of  { )  implies  in  this  case  the  strong 
convergence  of  the  discrete  stresses  (c^)  towards  (o)  in  (LP  (Q)),,X,,. 


a 


3  AUGMENTED  LAGRANGIANS 

3.1  Formulation  of  the  discrete  problems  as  saddle-point  problems. 

In  view  of  the  numerical  solution  of  the  approximate  viscous  flow  problem  (2.10)  by 
augmented  lagrangian  techniques/  we  must  first  reformulate  (2.10)  under  a  slightly 
different  form. 

lto  do  that,  observe  that,  if  we  replace  by  its  Ho(fl)  projection  over  the  space 

of  continuous  functions  whose  restriction  to  each  subelement  fi*  is  a  first  degree 
polynomial/  we  can  rewrite  (2.10)  as 


(3.1)  Minimize  F(D(w^))  +  G(x^)  over  1^,  with 


(3.2)  D(wh)  “V^V«h  +  V*?), 


f  :  Yh  *  R. 


(3.3) 


F(V  “  ^*<V  d*' 


\  * 


(3.4) 


G<V  “  -  -  Jr  9-h  da' 

1  2 


(3.5)  1^  -  {1^:0  ♦  RNXN,  *  V  E^i  ±  e  <Po(f2^ )  )N*N ,  V  i  -  1,  2N,  Vt  -  1,  . 

I0! 

If  we  follow  the  methodology  of  FORTIN-GLOW INSKI  [1982],  we  can  then  replace  (3.1)  by 
its  augmented  lagrangian  formulation 


(Find  a  saddle-point  ffvj  ,  V'  V  of  the  augmented  lagrangian 

WW  *  F(V  +  ®‘V  *  f,D(V  -  Vo, 2  -  <  V  D(V  -  s  > 

over  the  set  (K^  x  Y^)  * 
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where  R  is  any  positive  number  and  where  <• ,•>  denotes  the  classical  L  (ft)  scalar 
2  NxN 

product  over  (L  (fl))  .  Observe  that  (3.6)  imposes  the  Incompressibility  condition  on 

the  continuous  variable  (m^  must  belong  to  the  set  of  approximately 

incompressible  velocity  fields)  but  minimizes  the  nonlinear  functional  F(-)  with  respect 
to  the  piecewise  constant  variable  G^.  In  other  words,  there  is  a  splitting  of  the 
difficulties  of  our  problem  (nonlinearity  and  incompressibility)  between  these  two 
variables. 


THEOREM  3.1:  The  augmented  lagrangian  problem  (3.6)  and  the  approximate 


incompressible  viscous  flow  problem  (2.10)  are  equivalent:  to  any  solution  v^  o£ 

(2.10),  one  can  associate  a  solution  {{v^  ,  V  ,  Xh)  of  (3.6)  and  conversely.  Moreover , 

**h  Is  equal  to  D( )  and  there  exists  an  approximate  pressure  field  ph  ^n  such 

that  the  approximate  stress  tensor  field  ( -1^  -  p^  Id )  satisfies  the  discrete  equilibrium 

equations  and  constitutive  laws  (2.11). 

Proofs  Pirst,  let  v,_  be  a  solution  of  (2.10)  and  let  o.  and  p  be  the 
'  n  h  h 

associated  discrete  stress  and  pressure  fields.  From  (2.11),  Theorem  2.1,  we  have 


<ah  +  PhM>>Gh  «  Jn  {VW  -  W*  v  S  e  V 


j  o  *D(w  )dx  ”  ]  f ■  v,  dx  +  1  da,  V  *  e  V '  . 

q  h  h  a  h  r  h  hh 

2 


In  particular,  taking  as  (s^-v^)  where  ^  is  any  element  of  K^,  and  since,  by 

construction  of  Kj,,  div  (s^-^)  iB  ec!ual  to  zet°  m  the  dual  of  Ph,  we  have 


Jn  (tfhA  "’'"‘W*'  -  j  f*«VTh,dX  +  *  »*'VV  da'  V  *h  e  V 

2 


Substracting  this  to  the  first  inequality  of  (3.7)  yields 


wv  -  -  ?  ,D<v  -  vv2 .2 


for  any  {^,<^}  in 


x  V.  .  Since  in  addition 


'VW 


-I^Id)  -  yv*h 


’V 


J(V- 


V  U.  €  Y 


h' 


{{vh,^},-gh-|^M}  is  indeed  a  saddle-point  of  the  augmented  lagranglan  iR(*,*,») 
over  x  Yh>  x  Yh- 

Conversely)  let  { { v^ , } ,  X^)  be  a  solution  of  the  augmented  lagranglan  problem 
(3.6).  Then,  we  must  have 

VWV  >  WW’ 

which  can  only  hold  if  we  have 
(3.8)  ^  -  D(vh)  . 

Taking  (3.8)  into  account,  the  second  saddle  point  inequality  yields 


,3-9)  WD(V'V  *  VWV'  vfW  e  *h  *  V 

In  particular,  by  taking  Cj,  as  O(w^) ,  (3.9)  implies 

J(V  «  JtV  -  v\  e  V 

and  Vjj  is  indeed  a  solution  of  the  original  minimization  problem  (2.10). 

To  further  characterize  any  solution  ^vh,Hh^’  ^h^  of  the  augmented  lagranglan 
problem  (3.6),  we  again  use  (3.8)  and  (3.9).  From  (3.8),  is  necessarily  equal  to 
D(vh>.  On  the  other  hand,  introducing  the  space  defined  in  (2.13),  (3.9)  can  be 

rewritten  as 
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Equivalently,  if  we  consider  L^( v^+* , D(  >+• ,  X^)  as  a  convex  function  of  the  pair 
on  the  space  X^  x  Y^,  we  can  write  (3.9)  as 

{0,0)  e  3  !R(vh+0,D(vh)+0,Xh)  in  X^  x  Y*  . 

A  direct  calculation  characterizes  the  elements  of  this  subgradient  as  the  pairs 
{g^dj^}  of  x  Y^  such  that 


j  D(g)*D(»  )dx  <  -j  X*D(w)dx  -  j  f»w  dx  -  J  g* w  da ,  V»  e  X  , 
ft  ^  h  nhh  ft  h  r2*' 


j«  VV*  ‘  Jn  +  S’  -  +  VV*'  v  **,  e  V 


Setting  g^  and  to  zero  in  (3.10),  we  simply  obtain  the  variational  system  (1.9)- 

(1.11)  with  (a  )  =•  -X  .  As  seen  in  the  proof  of  Theorem  2.1,  this  in  turn  implies 

D  h  n 

(2.11)  with  oh  ■  -Xh  -  p^  Id,  and  our  proof  is  complete.  a 


RFWARK  3. 1 !  In  order  to  accelerate  the  convergence  of  the  algorithm  to  be  used  for 

the  solution  of  the  augmented  lagrangian  problem  (3.6),  it  is  usually  better  to  replace, 

2 

in  the  definition  of  the  augmented  lagrangian,  the  classical  L  (ft)  scalar  product  by  an 

equivalent  weighted  scalar  product  of  the  type 

<  C,  D  >  =  J  r(x)  C  •  D  dx. 
ft 

06 

Here  r(x)  is  a  strictly  positive  scalar  function  of  L  (ft),  bounded  away  from  zero, 
which  can  be  arbitrarily  chosen.  Proper  choices  for  this  function  will  be  discussed 
later.  With  this  new  scalar  product,  LR( •»•,•)  becomes 


WW  =  J  pi(V  *  ’  J  **  -  J  da 

1  2 

+  ^  j  r(x)  |<^-D(wh)|2  dx  -  J  r(x)uh*(D(«h)  -  G^dx  • 
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3.2  Numerical  algorithm  The  fundamental  interest  of  the  equivalent  augmented 
lagrangian  formulation  (3.6)  la  the  existence  of  a  very  cheap  and  simple  algorithm  for  its 
numerical  solution.  This  algorithm  combines  an  Uzawa  algorithm  for  the  solution  of  the 
saddle-point  problem  and  a  block-relaxation  technique  for  the  solution  of  the  minimi ration 
problems  associated  to  the  primal  variable  ,  (^ ) .  Dropping  the  subscript  h  from 
all  variables  for  simplicity,  this  algorithm  is 


(3.11)  Let  {X°,  H  be  given  in  x  Kh' 

Then,  for  n  >  0,  h”  '  and  Xn  being  known,  we  compute  {en,Hn}  ^n  x 
by  block-relation,  i.e.  by  setting 


and  by  computing  sequentially  vjj  and  by  solving 

<3-12’  <-i'A">  <  V*'<-i'xn)'  V 

(3.13)  lR(eR.  l£,  X")  <  LR(V^C,Xn),  VG8Yh, 

Once  {*",■"}  is  known ,  the  Lagrange  multiplier  X  is  updated  by 

(3.14)  X”  1  -  Xn  -  R(D(vn)  -  h"). 

Many  variants  exist  for  this  algorithm  and  are  described  for  example  in  FORTIN-GLOWINSKI 
(1982].  Usually,  the  block-relaxation  (i.e.  the  loop  ( 3. 12)-( 3. 13)  on  k)  is  only  carried 
out  for  one  to  five  iterations. 

Observe  that  the  above  algorithm  only  considers  one  variable  at  a  time  and  therefore 
takes  full  advantage  of  the  splitting  of  the  difficulties  achieved  by  the  saddle-point 
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formulation  (3.6).  First,  in  (3.12),  the  matrix  field  and  the  multiplier  ln  are 

supposed  to  be  known,  and  the  algorithm  minimizes  the  augmented  lagrangian  t R  with 
respect  to  the  velocity  field  w  in  As  function  of  the  velocity  field,  iR  is 

quadratic  and  corresponds  to  the  energy  dissipated  by  an  incompressible  Stokesian  fluid, 
flowing  viscously  under  the  action  of  the  external  loads  {f,g}.  In  other  words,  (3.12) 
is  a  classical  linear  stationary  Stokes  problem,  discretized  by  mixed  finite  element 
methods.  Many  numerical  techniques  are  available  for  its  solution,  and  we  refer  to 
TAYLOR-HOOD  (19731,  GIRACJLT-RAVIART  (1979]  or  GI/WINSKI-PIRONNEAU  [1979]  for  the  practical 
description  of  such  techniques.  In  our  numerical  experiments,  we  will  choose  a  conjugate 
gradient  method  operating  on  the  hydrostatic  pressure  space  P^,  which  only  requires  the 
inversion  of  sparse,  fixed,  positive  definite,  symmetric  matrices  and  therefore  only  uses 
little  computer  running  time  and  memory  core  (FORTIN-GLOWINSKI  [1982  p57]).  In  any  case, 
most  finite  element  codes  now  propose  efficient  subroutines  for  the  solution  of  the  Stokes 
problem,  which  can  be  blindly  used  for  solving  (3.12). 

Then,  the  algorithm  supposes  the  velocity  field  and  the  multiplier  ln  given, 

and  in  (3.13)  minimizes  LR  with  respect  to  the  matrix  field  G  in  Yh-  The 
incompressibility  condition  and  the  spatial  derivatives  of  G  are  not  involved  in  (3.13)i 
this  is  an  unconstrained  local  convex  minimization  problem  whose  numerical  solution, 
described  in  details  in  the  next  section,  reduces  to  the  solution  in  parallel  of 
independent  convex  minimization  problems  set  on  RN(N“2  or  3). 

Finally,  after  a  few  resolutions  of  (3.12)  and  (3.13),  the  algorithm  updates  the 
multiplier  Xn  by  the  explicit  formula  (3.14),  so  that  the  constraint  Dlv")  “  rf1  can  be 
better  satisfied  by  the  solution  of  (3. 12)-(3. 13) ,  and  then  returns  to  (3.12)  and  (3.13). 

3.3  Convergence  of  the  algorithm  ( 3. 1 1 )-( 3. 14 ) .  We  now  study  the  convergence 
properties  of  the  above  Uzawa  algorithm,  considering  the  basic  particular  case  where  only 
one  iteration  of  block-relaxation  is  done  at  each  step  of  the  Uzawa  algorithm.  In  our 
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-  ■  im 


at n? 


2 


study.  It  will  oe  most  f portent  to  worn  on  »h  witn  tne  precise  weiqnteq  L  (W)  norm 
which  is  used  in  ths  construction  of  the  augsMnted  lagrangian  LR  (sss  Remark  3.1). 

Then,  if  ws  dsnots  by  { v,  1,  X  }  ths  solution  of  ths  augmented  lagrangian  problem 
(3.6),  by  Xn  ths  multiplier  calculatsd  in  (3.14),  by  B"  ths  matrix  fisld 
calculatsd  in  (3.13)  snd  by  e"  ths  vsetor  fisld  calculatsd  in  (3.12),  we  can  prove 


CONVERGENCE  THEOREM  3.2i  Under  the  assumption  of  the  existence  Theorem  1.1  (01 
convex,  continuous,  coercive),  ths  sequence  {Xn}  is  bounded  in  Yh,  the  difference 
(Ddr")  -  rf1)  converges  to  taro  in  Y^,  and  the  quantity  Fl^1)  +  Gfv")  converges  towards 
the  dissipated  energy  rate  J(t).  If  in  addition  (2.16)  is  satisfied  together  with  the 
first  line  of  (2.17)  ( P1 ( • )  uniformly  convex  on  the  bounded  sets  of  Yh),  then  the 

sequence  {v  Srf1}  converges  towards  {e,B}  strongly  in  x  Yjj.  Finally,  if  the 
Internal  dissipation  potential  P1  ( • )  is  continuously  differentiable,  and  if  its  gradient 
is  invertible  with  a  coercive  and  Lipschlts  continuous  Inverse,  that  is  if  P^d) 
satisfies 


rSn  r(x)  fc,  -  C2|2dx  <  c23  /  r<«)  13^(8,)  -  3P,<ej)|2 


(3.15)  / 


J  (ap^c,)  -  ))•<«, -a2>am  >  c14  J  r<«) I -  aP^c^l2**. 


for  any  and  Gj  ijn  Yh,  then  we  can  prove  that  the  sequence  {»n,Bn,Xn}  converges 

linearly  towards  {v,H, X}  in^  x  Yh  x  Yh  with  an  asymptotic  constant  bounded  by 

2  V 

C15  -  ( 1-2RC14/(  1+C13R)  )  2  . 

Proof»  Since  {v,H,X}  is  a  solution  of  the  saddle-point  problem  (3.6),  the 
following  extrenality  relation  is  satisfied: 
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(3.16)  F(^)  +  G(t")  «■  <  X,  a"  -  IXv")  >  >  F(B)  ♦  G(»)- 

Moreover,  by  construction,  the  solutions  vn  and  H11  of  (3.12)  and  (3*13)  satisfy  the 

extreaality  relations 

G(w)  -  G(t")  +  <  R(D(Tn)  -  h"'1)  -  Xn.  !«•-▼")  >  >  0,  V  V  e  lt^, 

F(G)  -  RHn)  ♦  <  R  (Hn-D(Tn))  +  x",  H  -  Hn  >  >  0,  V  G  e  Y.  , 

n 

respectively.  By  addition,  setting  »  ■  »,  and  G  •  B,  we  get: 

(3.17)  F(B)  +  G(v)  -  RlD(v")  -  ^l2  -  R  <  Hn  1  -  b",  D(e-e")  > 

♦  <  X",  D( »" )  -  b"  >  >  F(Bn)  ♦  G(v")  . 

Adding  (3.17)  to  (3.16),  we  then  obtain 

(3.18)  -  RlD(wn)  -  Bn»2  -  R  <  b"-1  -  a",  Wt-f”)  >  ♦  <  x"-X  ,  IXv")  -  b"  >  »  0. 

Combining  (3.18)  with  the  construction  (3.14)  of  Xn+1  finally  yields 

(3.19)  lx"  -  XI2  -  IX"+1  -  XI2  >  R2  lD(wn)  -  b"|  2  +  2R2  <  b"”1  -  b",  IXv-w")  >  . 

On  the  other  hand,  using  (3.13)  and  (3.14)  at  iteration  (n-1),  we  can  estimate  the 
right-hand  side  of  (3.19)  by  standard  algebraic  manipulations.  Exactly  as  in  FORTIN- 
GLOWINSKI  (1982  p117,  equations  (5.17)  to  (5.24)1,  setting  PQ  -  0,  0  •  F  and  Inverting 

the  sign  of  X,  we  have  the  following  estimate 

2R2  <  D(vn-e),  Bn  -  Bn‘1  >  >  R2(IBn-BI2-IBn'1-BI2)  +  R^h"  -  h""1!2, 

which,  combined  with  (3.19),  gives 


-28- 


(3.20)  <IXn  -  XI2+R2IHn"1-BI2)  -  (IXn'M-X|2+R2IBn-BI2)  > 

R2I  D(v")  -  rfV  ♦  R2 1  rf1  -  b""1!2  . 

The  positive  sequence  IX  -XI  ♦R2IBn  -Bl2  ia  therefore  decreasing  and  thus 
converges  to  a  Unit.  This  implies  that  the  right-hand  side  of  (3.20)  nust  converge  to 
sero  and  wa  finally  obtain 


IXn  -  XI2  +  R2IB°  ^  -  Bl 2  is  bounded, 

lin  I D( v" )  -  rf4!2  -  0, 
n+t" 

lin  IB14  -  b"*1!  »  0  . 

TYi«ae  convergence  results,  used  back  In  (3e16)  and  (3*17)  obviously  Imply  the  convergence 
of  F(rfM  +  aV1)  towards  F(B)  +  G<v). 

Now,  if  01  ( • )  ia  uniformly  convex  on  the  bounded  sets  of  Yh,  the  convergence  of 
the  energy  rate  and  the  boundedness  of  B"  imply  the  convergence  of  the  arguments  tf1 
and  v*  respectively  towards  B  and  v. 

Finally,  since,  from  the  Korn's  inequality,  BTB  is  an  isomorphism  from  onto  its 

dual,  we  can  prove  the  linear  convergence  of  the  sequence  {v™,  b",  Xn}  by  applying  a 
result  of  LIONS -MERC  IER  (1979,  Prop.  4,  p  970]  which  will  be  applicable  here  as  soon  as 
(3.15)  is  aatisfied(  see  FORTIN-GLCWINSKI  [1982,  p  300]  for  more  details).  q 

R  PI  ARK  3,1i  Condition  (3.15)  is  satisfied  at  least  locally  for  Norton  materials.  It 
is  not  satisfied  in  the  general  case,  but  linear  convergence  of  the  sequence  {v11,  b") 
can  still  be  observed  numerically  in  almost  any  case.  q 

REMARK  3.2i  The  asymptotic  constant  appears  numerically  not  to  be  optimal. 

Nevertheless,  its  expression  as  a  function  of  C13  and  C14  will  indicate  the  right 
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1 


strategy  to  follow  for  the  choice  of  the  parameter  R  and  of  the  weight  r(x).  Since 
Cjj  and  highly  depend  on  the  weight  r(x)  used  in  the  definition  of  the  augmented 

lagrangian  I.R,  i«  a  function  of  R  and  of  r.  TTie  right  strateg/  for  the  choice 

of  R  and  of  r  now  consists  in  trying  to  keep  C15  as  small  as  poss  ble.  By 
choosing  R  close  to  1/C13,  becomes  approximately  equal  to  »  (  1  -C 14/2c ^ j ) 

By  taking  the  weight  r(x)  so  that  the  products 

(3.21)  J  OP  (■  )  -  3P  ,(■-))  •  <C  -G  >dx  and  (3.22)  J  r(  x)  (6,-C,  )•  (■,  -H,  )dx 
fl  i  1212 

remain  close  to  each  other  when  and  are  in  the  neighborhood  of  H,  the  ratio 

C14^C13  9etB  close  to  1  and  reaches  the  value  1//T  .  The  final  strategy  for 

choosing  R  and  r  is  therefore: 

(1)  choose  r(x)  by  matching  (3.21)  and  (3.22)> 

(ii)  take  R  close  to  1/Cjj,  which  will  usually  be  close  to  1  if  r(»)  is 
properly  chosen. 


REMARK  3.3:  If  P1  is  quadratic  and  if  we  have  equality  between  (3.21)  and  (3.22), 
then  for  R  -  1,  w”  converges  in  2  iterations  and  rf1  converges  linearly  with 
asymptotic  constant  .5  ( FORTIN-GLOWINSKI  (1982,  p  119)).  In  other  cases,  with  proper 
choices  of  R  and  of  r,  we  usually  observe  linear  convergence  of  ( »" , }  with  an 
asymptotic  constant  around  .7.  O 


REMARK  3.4:  Linear  convergence  compares  unfavorably  to  the  quadratic  convergence 
expected  for  conjugate  gradient  or  for  Newton  algorithms.  But  Newton  method  requires 
Vy  to  be  twice  differentiable,  the  factorization  of  quite  a  few  finite  element  matrices, 
and  its  convergence  rate  can  be  very  slow  for  weakly  convex  dissipation  potentials.  In 
general,  it  is  not  a  good  method  for  solving  (2.10).  Oi  the  other  hand,  a  conjugate 
gradient  method  with  preconditioning,  of  the  type 
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1 


I 


•  taka  in  K^i 

•  aolva  <0(9),  D(w)  >*<J'(m),  w>,VweK» 

o  on 

•  «et  s  *  9  / 

o  o 

•  for  n  -  0,  until  satisfied  do 

p  ■  Arg  Min  J(u  -p«  ), 
n  n  n 

a  x  ■  n  -  p  s  , 

n+1  n  n  n 

solve  <  Dtg^),  D{«)  >  -  <  J'(«n+1),w  >  ,  V  v  e  1^, 

Yn  '  <  “Vi1,  WVA>  >  7  <  °<V'  «n  >  ' 

Vi  “  Vi  *  Tn  *n' 
end  loop  on  n  ; 

2 

where  <•,  •>  la  an  adequate  weighted  L  (fl)  scalar  product  on  Yh,  will  only  be 
efficient  if  t^(*)  is  differentiable,  if  the  scalar  product  on  is  correctly  chosen 

and  if  a  very  efficient  Stokes  solver  is  available  for  coeiputing  gn+1>  If  this  is  the 
case,  the  conjugate  gradient  method  will  be  twice  as  fast  as  the  Uzawa  algorithm  (3.11)— 
(3.14).  If  this  is  not  the  case.  Algorithm  (3. 1 1) -< 3. 14)  appears  to  be  one  of  the  only 
reasonable  numerical  method  for  solving  (2.10). 
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«.  THE  PROBLgtS  IN  DEFORMATION  RATES. 

4.1  The  local  problems.  Problem  (3.13)  appears  as  one  step  of  the  algorithm 
proposed  herein  for  the  numerical  solution  of  the  viscous  flow  problems  In  quasistatic 
viscoplasticity,  once  these  problems  have  been  approximated  by  simplicial  finite  elements 
of  order  one  and  decomposed  under  an  au^nented  lagrangian  form.  Recall  that  here,  this 
problem  consists  in 

(3.13)  Minimizing  |.r(t,*,X)  over  Yh 

with 

L  (v,G,X)  -  J  P?(G)  dx  -  J  f*v  dx  -  J  g*v  da 

r  n  n  r2 

•f  2  J  r(x)|c-(Vv+VvT)/2|2dx  -  ]  r(x)X*  [  (Vw+VvT>/2  -  g]  dx. 

2  n  Si 

Yh  -  {O:  fl  ♦  ,  D,  i  e  <Po<nJ )  J1****,  Vi  -  1,  2N,  vt  -  1,  sh), 

'nt 

and  that  it  is  the  only  nonstandard  step  in  this  algorithm,  the  other  steps  consisting  of 
linear  Stokes  problems  and  explicit  variables  updating,  respectively. 

Since  all  the  elements  of  Yh  are  matrix  fields  which  are  constant  on  each  8^,  and 
since  the  functional  L„  does  not  involve  any  distributional  derivative  of  G,  Problem 
(3.13)  can  equivalently  be  written  as 

(4.1)  V  i  -  1,  2N,  V  t  -  1,  1^,  Minimize  J*(G)  over 
with 

(4.2)  J^(6>  -  C?(G)  +  £2  |G| 2  -  rG  •  <R(?v*VvT>/2  -  X).  . 

1  '  2  In; 

Here,  we  are  simply  using  the  fact  that  the  minimum  value  of  the  sum  of  independent  terms 
is  equal  to  the  sum  of  the  minimum  value  of  each  term.  Then,  Problem  (3.13)  reduces  to 


A 
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II 

th«  solution  in  parallel  of  N^x2  local  indapandant  convax  minimi ration  problaaa  act 

on  rNxN  (N  -  2  or  3). 

a 

Using  a  general  purpose  minimisation  algorithm  for  the  solution  of  each  local  problem 
(4.1)  is  not  adviseable  here  for  two  main  reasons: 

(i)  such  an  algorithm  is  very  difficult  to  implement  because  it  must  be  able  to 
handle  general  nondif ferentiable  convex  dissipation  potentials 

(ii)  such  an  algorithm  is  usually  expensive  in  computer  running  time. 

An  easier  and  more  efficient  strategy  consists  in  adapting  each  time  the  minimisation 
algorithm  to  the  specific  class  of  potentials  V 1  which  is  under  consideration.  Doing 
that,  we  have  most  of  the  times  been  able  to  reduce  each  local  problem  (4.1)  to  a  one¬ 
dimensional  convex  minimisation  problem  set  on  The  remainder  of  this  report  will 

describe  the  derivation  of  such  efficient  numerical  techniques  in  the  case  of  Norton 
materials,  of  Bingham  materials,  and  of  Tresca  type  materials  in  plane  stresses, 
respectively.  But  before,  we  will  derive  a  very  useful  simplification  of  the  local 
problem  (4.1). 

4.2  Reduction  of  the  local  problems.  He  begin  by  recalling  several  well-known 
results  of  matrix  theory,  which  will  enable  us  to  reduce  the  local  problems  (4.1) 
which  are  set  on  to  local  convex  minimi  gat  Ion  problem*  set  on  R**,  (N  ■  2  or  3)  • 

Imams  4.1  (VON  NEUMANN!  1937) ) .  Let  A  and  B  be  two  matrices  in  *NxN  with 
singular  values  >  •••  >  >  0  and  8^  >  B2  >  •••  *  8^  >  0.  Then 

N 

Tr(A  B)  <  l  a  8  . 

i-1  1  1 
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Lemma  4.2.  Let  A  and  B  be  two  symmetric  matrices  In  RNjjN  with  eigenvalues 


A  >A_>***>A  and  B.  >  B  A  •••  >  B  .  Then 
1  2  N  -  1  2  N  - 


A  •  B  -  Tr  {A  B)  <  l  A.  B 
1-1  X  1 


Proof :  The  result  follows  from  the  decomposition 
Tr  (A  B)  -  Tr[(A-AN  Id)  [B-Bn  Id))  +  A^Tr ( B)  +  B^Tr  ( A)  -  N  A^ 
and  from  the  application  of  Lemma  4-1  to  the  first  term  of  the  right-hand  side.  ^ 


Lemma  4.3  Let  D°  be  a  diagonal  matrix  with  diagonal  terms  >  •  ■ 

and  A  be  a  symmetric  matrix  with  eigenvalues  A ^  >  A^  >  • • “  >  A^ .  Then 


>  D. 


N 

Max  [Tr(PTDd  P  A)  j  -  Tr  [0  Dd  g1*)  “  l  DA, 

T  i-1 

PP  -Id 


where  Q  is  an  orthogonal  matrix  which  diagonalizes  A  with 

(QT*0,ii  -  a±  . 

Proof :  For  p  given,  let  B  be  the  matrix  defined  by 

T  d 

B  -  P  D  P 

and  whose  eigenvalues  are  D, .  TTien ,  from  Lemma  4.2,  we  have: 


_T_d_ 


N 


Tr  [P  D  P  A)  -  Tr  (AB)  <  l  A  D  . 

1-1  1  1 

On  the  other  hand,  we  also  have 


N 

Tr[QDdOTA]  -  Tr  (gTQDdgTAg)  -  Tr  [DdgTAg)  -  l  A  D 

i-1 


and  the  result  follows. 


We  are  now  ready  to  prove  the  main  result  of  this  section,  which  reduces  the  local 
problems  to  convex  minimization  problems  set  on  R^(N  -  2  or  3). 
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THEOREM 


4.  1 :  It  the  internal  dissipation  potential  ( • )  is  convex  and  isotropic. 


then  the  solution  H  of  the  local  problems  (4.1)  is  given  by 


(4.3) 


matrix  Aj  (the  first  column  corresponding  to  the  biggest  eigenvalue  and  so  on)  and 
where  Bd  is  the  diagonal  matrix  of  ttt,xN  solution  of 


(4.4)  Minimise  Ip^D1)  +  I  Dd  | 2  -  )  Ai<D<1'ii^  on  °N 

i»1 


T  d 

G  =  P  D  P  , 


Above  DN  denotes  the  space  of  diagonal  matrices  of  RNxN  and  ( A^ )  denotes  the  set  of 
eigenvalues  of  A*  (A;  >  A^  >  •••  >  A^) . 

NxN 

Proof.  First  observe  that  any  matrix  G  of  R  can  be  decomposed  into 

(4.5) 

where  P  and  Dd  are  two  independent  matrices  of  R*IxN,  P  being  orthogonal  (P  PT  * 

Id)  and  being  the  diagonal  matrix  whose  diagonal  elements  are  the  eigenvalues  of 

N  N 

G.  Then,  if  we  denote  by  0  (respectively  D  )  the  space  of  orthogonal  (respectively 
diagonal)  matrices  of  RNx  ,  the  local  problem  (4.1)  becomes 


(4.6)  V  i,  Vi,  Minimize  Jd(P,Dd)  over  0N  x  1)N, 


with 


Jj(P,Dd)  -  J*(G)  =  P*(PTDdp)  +  I^IpV1?!2  -  (PTDdPl  •  A*. 


But  since  V^(‘)  and  |*|  are  isotropic  (P*(P  D  P)  =  P*(D  )),  we  can  rewrite  J*  as 


J^(P,Dd)  =  P®(Dd)  +  y-  |Dd|2  -  (PTDdP)  •  Aj. 
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Now,  let  be  the  solution  of  (4.4)  (unique  since  the  function  to  minimize  is 

strictly  convex)  and  let  PH  be  the  permutation  matrix  which  reorders  the  diagonal 

d  „  e 

elements  of  in  the  decreasing  order.  Since  A  ^  >  •  •  •  >  and  s*nce  ”  ^  (  *  ) 

isotropic,  we  have 


r»e,J Pjd_  .  rR 

V  P  H.  P  )  +  — 
V  H  £  H  2 


N 

-  1 

i=  1 


T  d 

VWh'i  ^ 


_d,2 


N 

-  1 

i=  1 


A.  (n£)  . 
1  l  1 


But  Hd  is  the  only  solution  of  (4.4),  thus  pT  Hd  PH  has  to  be  equal  to  H^,  which 
£  H  £  £ 

d 

means  that  the  diagonal  elements  of  are  already  placed  in  a  decreasing  order.  From 

Lemma  4.3,  this  implies 


(4-7)  ^  VHi>i  '  H*  °I>  *  \  * 


Moreover,  since  H  is  solution  of  (4.4),  we  also  have 
l 


P*(Hd)  +  ~|  |  -  l  A.(H^).  <  0*(Da!  +  f-  1 0°  |  ^  -  l  A^.,  V  D°  e  d", 

i«1  i=1 


e  ,  d  Rr  •  d  i  2 


and  in  particular,  if  Pp  is  the  permutation  matrix  reordering  Da 


p>?>  ♦  ^  '-d 


e , _T  d 


i-r  '  -  'V  -  }  VV i  4  p,(pddV  +F 

1=  i 


*  'pyp0|2  -  (pVpD)iAi.  V  Dd. 


Now,  from  (4.7),  Lemma  4.3,  and  the  isotropy  of  p^*),  this  implies 


+F  l"tl2  -  <wl>  •  <**> 


e  d  Rr  *  d  1 2  ,  T  d 

<  PJD  >  ♦  1°  I  -  ( P  D  P) 


(A,) 


for  any  Dd  in  PN  and  any  P  in  0N.  So,  finally,  we  get 


i  d  i  d  f  d ,  N  N 

J£  (0£#Hf)  *  *VP'D  >  vfp'D  )  .  in  0  x  P 


is 
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,  d, 

In  other  words,  (Q^,  Hf I  Is  a  solution  of  (4.6).  By  construction,  this  implies 
that  H,  given  by  (4.3),  is  a  solution  of  the  original  local  problems  (4.1).  Our  proof  is 
therefore  complete  because,  since  J3  is  strictly  convex,  such  a  solution  is  unique.  q 

REMARK  4. 1 .  It  is  well  Known  that  the  internal  dissipation  potential  is  isotropic 
and  convex  for  all  standard  isotropic  viscoplastic  solids  and  all  standard  viscoplastic 
fluids.  Theorem  4.1  can  therefore  be  applied  in  most  practical  situations.  ^ 

RIMARX  4.2i  In  Theorem  4.1,  the  relation  (4.3)  simply  expresses  that  Hi  and  A* 

K 

have  the  same  eigenvectors.  In  essence,  this  is  the  discrete  equivalent  of  the  well  known 
result  which  states  that  principal  stresses  and  principal  strains  are  parallel  in 
isotropic  elasticity  or  visoplasticity .  a 


S.  NORTON  VISCOPLASTICITY 

5.1  The  local  problems.  Norton  viscoplasticity  corresponds  to  one  of  the  easiest 
possible  case  where  the  Internal  dissipation  potential  P  (•)  is  given  by 

(5.1)  0(G)  -  -  (k/2)P|G|P  -  -  (k/2)P  (  l  G*  !P/2  . 

P  P  i,j-1  3 

Each  local  minimization  problem  (4.1)  now  becomes 


(5.2)  Minimize  J*(G)  over  RN*N  with 

— t  s  — - 

r  l6)  “  pW)«)P  +  f-  |g|2  -  A*  •  G, 

< 

-  rf  |  (^v+VvT)  -  l)i  ,  , 

L  K 
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whose  solution  is  given  by 


THEOREM  5.1:  The  solution 


of  the  local  minimization  problem  (5.2)  is  of  the 


form 


(5.3) 


i  i  ii 


where  x  is  the  solution  of  the  one-dimensional  convex  minimization  problem 


15.4)  Minimize  {-  (k/2)P  ^  ~  y2  -  I I y }  over  R.  . 

p  2  x  ▼ 


,  i ,  i  i  NxN 

Proof :  Suppose  that  the  norm  | |  of  the  minimizer  of  over  R  g  is 

known.  Then,  H*  minimizing  J*  over  R*^N  will  in  particular  minimize  J*(-)  over 

the  set  of  matrices  of  with  fixed  norm  | | .  But  this  last  minimization  problem 

reduces  to  the  maximization  of  the  scalar  product  A*.  G  over  a  Bphere  of  RN*N  and  ^ts 


solution  is  given  by 

(5.3) 


hJ  -  aJ  !<)/)<!* 


By  plugging  (5.3)  into  the  expression  of  J^,  the  minimization  of  over 

finally  reduces  to  finding  the  unknown  norm  | (  which  has  to  mlninlmlze  J^lA^y/lA^I) 
over  the  set  R+  of  positive  numbers.  This  last  problem  is  precisely  (5.4)  and  our  proof 
is  complete .  q 

In  practice,  the  numerical  solution  of  each  local  problem  (5.2)  uses  Theorem  5. 1  and 
is  achieved  by; 


a)  the  computation  of  the  solution  | |  of  the  one  dimensional  convex  minimization 
problem  (5.4) 


b)  the  computation  of  by  the  explicit  formula  (5.3 


.3)  . 
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i 

f 


The  numerical  solution  of  (5.4)  can  be  achieved  for  example  by  using  the  one  dimensional 
Newton  algorithm  below 


data:x°  ■  solution 

of  (4.4)  at  previous  iteration; 

initialisation 

X  ♦ 

o 
x  ; 

j  «• 

1; 

repeat  : 

j  -  j 

+  1; 

<J  ♦ 

(k/2)P  xP_1  +  rRx  -  +  A*|/2; 

test 

:  if  |g|  below  tolerance  exit; 

dg  ♦ 

<p-1)(k/2)P  xP~2  ♦  Rr ; 

X  ♦ 

-30  -1 

max  (10  ,  x  -  (dg)  g); 

exit 

KJ ' 

x. 

Putting  together  all  the  steps  which  permit  the  numerical  solution  of  the  decomposed 
approximated  viscous  flow  problem  (3.6)  by  the  algorithm  ( 3. 1 1 )-( 3. 14 ) ,  we  finally  obtain 
the  simple  and  easy  to  code  computer  flow  chart  of  Figure  5.1. 


5.2  Numerical  results.  The  first  numerical  test  considers  a  horizontal  cylindrical 


hose,  with  external  radius  1.,  which  is  glued  on  a  rigid  core  on  its  internal  face,  and 
which  is  subjected  to  its  own  weight,  this  situation  may  represent  for  example  the 
cooling  process  of  the  plastic  coating  of  an  electrical  wire,  for  which  manufacturers  must 
verify  that  the  deformations  undergone  by  the  coating  during  cooling  remain  small.  In  a 
first  approximation,  strains  are  assumed  to  remain  plane,  and  the  coating  is  supposed  to 
be  made  of  an  homogeneous  Norton  material  with  k  -  .47,  p  «  1.4,  and  volumic  weight  .1. 

For  symmetry  reasons,  only  the  right  half  of  the  section  is  considered)  225  nodes  are 
used  to  approximate  the  velocity,  65  nodes  are  used  for  the  pressure.  Only  one  block- 
relaxation  iteration  is  done  at  each  Uzawa  step  and  after  60  iterations  of  the  Draws 
algorithm,  the  error  Ih"  -  D(v”)l  has  decreased  by  a  factor  of  10  7  (in  fact,  30 
iterations  were  more  than  sufficient  to  obtain  a  very  accurate  velocity  field) .  the  total 
CPU  time  was  approximatively  3mn  on  the  VAX  760.  Figure  5.2  represents  the  computed 
velocity  field  (magnified  7  times).  Ttie  shape  of  the  hose  after  one  second  of  flow  (the 

deformations  being  multiplied  by  40)  is  Indicated  on  Figure  5.3,  together  with  the  mesh 

used  for  the  pressure. 

In  the  entire  computation,  the  parameter  P  was  equal  to  1.  The  weight  r(x)  was 

i  20 . p-2 

equal  to  1  during  the  first  20  iterations,  then  updated  by  the  formula  rl«)  «  |K  I 

(see  Remark  3.2  for  justification),  kept  that  way  until  iteration  40  where  it  was 

i  40 i p-2 

finally  updated  by  r(x)  -  ]H  I  . 


6.  BINGHAM  VISCOPLASTICITY. 

6.1  The  local  problems.  Bingham  viscoplasticity  corresponds  to  an  internal 
dissipation  potential  ( * )  of  the  type 
(6.1)  P^G)  -  u|c|2  +  g/2  |Gi  . 

Each  local  minimization  problem  (4.1)  now  becomes 


RAPPORT-  OlOOO 
bnosam  rum  »  A  CATOT 
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(6.2)  Minimize  f (|S  +  M) |g| 2  +  q/2  |c|  -  A^  •  g)  over  RNxN 

2  l  - -  3 

whose  solution  is  simply  given  by  the  explicit  formula 

(6.3)  -Max  {0, 1-/5  g/|A*|}  aJ  /<rR+2u)  . 

In  practice,  the  whole  program  solving  the  flow  problem  (3.6)  for  Bingham  fluids 
still  corresponds  to  the  computer  flow  chart  of  Fig.  5. 1,  each  local  problem  being  now 
solved  by  (6.3). 

6.2  Numerical  result.  We  consider  herein  a  Bingham  fluid  flowing  viscously  through 
a  cavity.  Fluid  enters  at  the  upper  right  of  the  cavity  and  exits  at  the  upper  left,  with 
an  imposed  velocity  of  3.0.  No  slip  boundary  conditions  are  imposed  elsewhere. 

Dimensions  of  the  cavity  are  1.  for  the  main  square  and  .1  for  the  entrance  and  exit 
tubes.  The  fluid  viscosity  u  is  assumed  to  be  .01  and  velocities  are  supposed  to  remain 
plane . 

The  finite  element  mesh  uses  419  nodes  for  velocities  and  166  nodes  for  pressures. 

One  block-relaxation  is  done  at  each  Uzawa  step  and  we  take  R  -  .1  and  r(x)  -  1.0. 
Figures  6.1  and  6.2  represent  velocity  obtained  after  40  Iterations,  the  plasticity 
threshold  g^2  being  respectively  of  u  and  of  lOy.  As  expected,  the  domain  where  the 
fluid  is  at  rest  is  bigger  in  the  latter  case.  The  computation  time  for  each  case  was 
approximatively  10  mn  on  the  VAX  780. 

7.  TRESCA  TYPE  VISCOPLASTICITY  IN  PLANE  STRESSES. 

7,1  The  local  problems.  In  plane  stresses,  the  body  under  consideration  is  supposed 
to  be  very  thin  along  Xj  and  is  loaded  in  its  plane  so  that,  in  a  first  approximation, 
all  stresses  along  x^  are  equal  to  zero.  It  is  then  possible  to  eliminate  the  Xj 
direction  and  to  reduce  our  original  problem  to  a  two-dimensional  one  whose  domain  will  be 
the  middle  plane  section  of  the  body  and  whose  unknowns  will  be  the  in-plane  velocities. 
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These  inplane  velocities  need  no  longer  be  incompressible  since  any  reduction  of  the  plane 

section  can  be  compensated  by  a  corresponding  thickening  of  the  body.  Therefore,  in 

Sections  2  and  3,  K.  is  everywhere  replaced  by  v.  +  u  and  in  particular  (3.12) 
n  ho 

becomes 


Minimize 


l_(*.G,X)  over  V.  +  u  , 
R  n  o 


■ 

which  is  a  classical  linear  elasticity  problem  with  a  zero  first  Lame  coefficient.  As  for 
(3.13),  its  formulation  is  unchanged  and  it  still  reduces  to  the  local  problem  (4.1). 

In  plane  stresses,  a  Tresca  type  viscoplastic  material  corresponds  to  the  internal 
dissipation  potential 

(7.1)  pf(G)  <k/5)P  {Max(  |  | ,  |  | ,  |  G^+G^  | )  }P 

where  and  G2  are  the  eigenvalues  of  the  2x2  symmetric  matrix  G.  This  potential 

is  a  symmetric  convex  function  of  the  eigenvalues  of  G,  therefore  is  isotropic  and  convex 
(HILL( 1970] ) .  Moreover  it  satisfies  the  inequalities  (1.6)  for  coerciveness  and 
continuity.  On  the  other  hand,  this  potential  is  not  Btrictly  convex  and  not 
differentiable,  which  clearly  appears  in  Fig.  7. 1,  where  the  level  lines  of  0,  are 
drawn. 


The  viscous  flow  problem  associated  to  this  potential  (7.1)  can  still  be  solved  by 
the  Uzawa  algorithm  (3.11)— (3. 14 ) »  (3.12)  is  a  linear  elasticity  problem  and  (3.13) 
reduces  to  local  problems  whose  solution  H  is  given  by  (see  Theorem  4.1). 


(7.2) 


At  *  r,f  'V^V*T>  -  x>jni  •  -  i0  h2 


.  H  o 
H?  >  f  1  ..  )  , 


<  2 
'  j  (H1,H2)  <  j  wvv2).  v  e  F  , 


j  (D1,D2)  "  p  <k^)P[Max  {|di!,|d2!,|oi+D2|}|P  +  f2  (f,  +  -  a,  d,  -  a2  d2 
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Above,  as  before 


is  an 


Consequently,  the  computer  flow  chart  associated  to  the  Uzawa  algorithm  ( 3 . 1 1 )- ( 3 . 14) 
for  the  numerical  solution  of  viscous  flow  problems  in  Tresca  type  viscoplasticity  in 
plane  stresses  is  the  one  described  in  Figure  5.1,  but  with  the  solution  of  the  local 
problems  being  achieved  by  the  four  steps  above.  Among  these  steps,  only  one.  Step  3,  is 
not  explicit  and  its  solution  is  described  in  the  next  paragraph. 


7.2  Solution  of  Step  3.  By  definition  of  the  subgradient,  Step  3  is  equivalent  to 
(7.3)  {0,0}  e  9j(H1,H2). 

Therefore,  to  solve  Step  3,  we  first  begin  by  computing  the  subgradient  of  ;)(•,•)  over 
Since  we  know  from  Theorem  4.1  that,  at  the  solution  of  (7.3),  is  greater  or 

equal  to  Hj,  we  restrict  ourselves  to  the  half  plane  TTten,  a  direct 

calculation  gives: 

Case  1i  if  {D j ,D2 }  e  -  {{x,y>  e  RZ,  x  >  y,  y  >  0),  then; 


j(D1,D2)  "  p  (k/2>P(D1+D2)P  *  F  <D1+D2)  ‘  A1D1  *  A2D2' 

3  j  f  D1  ,D2>  -  {(u^Uj),  uA  -  (k/IjPfD^D^P'1  +  rR  DA  -  A^}  . 

Case  2 1  if  {D^.D,,}  e  K2  “  e  r2<  *  *  °«  Y  ”  <  then  < 

j(Dl'D2)  “  ^<k/2>P<D1)P  +  Di  '  A,0,' 


SjlD^Dj)  -  f(u1(u2),  u,«  k/ilD^P  +  rR  Dt  -  Ar-  A2  <  Uj  <  (k/2)P(D1)P'  -  A2) 
Case  3 1  if  {Dj.Dj}  e  K3  ”  Hx,y}  e  *2 ,  x  >  -y,  y  <  0},  then 

3<D1'D2)  ’  p  (k/Z  °1)P  +  I2  (01+D2!  ’  A1D1  '  A2D2’ 

3j(D1tD2)  -  {(u1,u2),  u1  -  (k/IjPtD^P-1  ♦  rR  D1  -  A^  . 

Case  4;  if  fDi'D2-'  e  K4  “  {fx,y}  e  *2'  x  “  K'  Y  <  °3'  then 

■j(D1,D2)  -  j  <-k/2D2)P  ♦  f  (DZ+DZ)  -  A,Dl  -  A2D2, 

3j(D,,D2)  -  {(u1(u2),  rRD,-  A,  <  u,  <  ( k/2 )P(D  1  )P_  1  +  rRD,  -  A,, 


rRD2  '  R2  ■  <k/2)P(-D2)P_1  <  u2  <  rRD2  -  A.,} 


Case  5:  if  ^di'd2^  e  *5  =  f (x,y)  e  R  ,  -y  >  x  >  0},  then 

f  j(Dl(D2)  =  1  (k/F)P(-D2)P  +  f  (d]+D22)  -  A,D,  -  A2D2, 

3j<t>  <D  >  »  u.  =  -(k/2)P(-D  )P  1  5  +  rR  D  -  A.). 

12  1  2  1  2  12  ii 

Case  6:  if  {D.,D„}  e  K,  =  {{x,y}  e  R2,  x  =  0,  y  <  0),  then 

-  -  l  Z  b  - 

j(D  ,D  )  =  -(k/2)P(-D  )P  +  *“  D2  -  AD, 

12p  2  22  2  2 

3j(D1,D2)  -  {(u1,U2),-(k*'2)P(-D2)P_1-  A^  -A^u^  -(k/2)P (-D2 )P~ ’+  rRD2“A2)  . 

Case  7:  if  {D^ ,D  )  e  =  {fx,y}  e  R2 ,  0  >  x  >  y) ,  then 

j<Dl'D2’  =  •;<I'/2)P(-D1-D2)P  t  (D2tD2)  -  A1D1  -  A2D2, 

3j(D1,C2)  *  {(Ul(u2),  UjL  -  -(k/2)P(-D1-D2)P_1  *  tR  D,  -  Ai>. 

Observe  that  the  sets  form  a  partition  of  the  half  plane  D1  >  D2>  Then,  by 

definition  of  (7.3)  and  since  its  solution  {H  ^ ,  H2 }  belongs  to  this  halfplane  (Theorem 
4.1),  we  have  : 

7 

(7.4)  (H-|  ,H2 )  =  O  { 1  e  Ki'  {0'05  e 

Therefore  the  solution  of  (7.3)  is  simply  the  solution  of  one  of  the  local  subproblems 
(the  one  which  admits  a  solution  for  the  given  data  of  A^  and  Aj) 

{0,0}  e  3j(xi,yi),  (xi,yi)  e  . 

Then,  once  that  for  each  subproblem  the  conditions  which  guarantee  the  existence  of 
solutions  are  explicited  and  that  the  algebraic  expressions  of  these  solutions  are 
computed,  is  simply  obtained  by: 

(i)  finding  which  subproblem  has  a  solution  for  the  given  values  of  A1  and  A2, 
checking  successively  the  admissibility  requirements  (conditions  for  existence  of 


by 


solutions)  of  each  subproblem: 


<  1 1 )  setting  e<3ual  to  the  corresponding  solution. 

Here,  these  computations  are  easy  to  carry  out  since  we  have  just  computed  the  algebraic 
expressions  of  <*!(•,•)  on  each  subset  For  example,  for  i  -  1,  we  have 

local  subproblem: 

(k/21  (x+y)P_1  +  rRx  -  A1  »  0, 

)  (k/2)P(x+y)P*1  +  rRy  -  A2  -  0, 
x  >  y  >  0; 


admissibility  requirement  (necessary  and  sufficient  condition  for  existence  of 
solutions ) : 


»2  >  t(A1  -  AjJ/rR] 


solution 


,x1  -  (z+fA^A  )/rR)/2,  yl  -  (z-(A1-A2)/rR)/2, 


z  minimizes  f—  (k/2)PtP  *  t2  -  (A.+A.Jt)  over  R  . 

— p  2  i  2  1  + 
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All  computations  done,  the  solution  of  Step  3  Is  finally  given 


for  A2  >  ((A^AjJ/rR)51-1  < 


H ?  -  U+fA^A^/rRl/a, 


H2  =■  [z-(ArV/rR)/2' 


z  minimizes  {— (k/2)PtP+  ~t 2-(A+A_ ) t }  over  R  , 

p  2  12  — — + 


for  <(A  -A„)/rR) 


for  A1  >  (-A2+(-A2/Rr)P_1)  >  0 


A  >  0  < 


H.  minimizes  {—(k/2)PtP+  r— t2-  A.t}  over  R 
*  P  2  i  + 


«2  -  0  I 

H.  minimizes  {—(k/2)PtP+  r^t2-A_t)  over  R 
l  —  p  2  I  + 


V. 


H2  -  Aj/Rr; 


for  <-A2+(-A2/rR)P_1)  >  Af  >  0 


(h  minimizes {-(k/2)PtP+  Rrt2  -(A  -A  ) t> 


1  «2,v,  over  R+ 


and  (A^rR)1*-1  +  Aj  >  -A2 


< 


vH2  "  -H1» 


for  -A2  >  A1  +  (A  ^/rH) P_1  >  0 


H1  -  A^/rR 


f  — H_ )  minimizes{~(k/2)PtP+Rrt2+A_t}  over  R  i 

<L  '  P  2  — —  + 


for  0  >  At  >  -( (A1-A2)/rR) 


p-1 


H1  '  ° 


( -H_ >  minimizes  {— (k/2)PtP+?^t2+A_t1  over  R 

2  "  p  2  2  ♦ 


f  H,  «  (z  +  (A1-A2)/rR)/2, 


for  -((ArA2,/rR,p-’  >  A,  <  H2  "  <-<VV/rR,/2' 


i  {  — z )  minimizes  {-( k/2)ptP-t~t2+<  A  ,+A.  )t]  over  R  . 
V  1  '  -  p  2  12  ♦ 
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In  the  above  formulae,  the  minimization  over  R+  is  numerically  achieved  by  using  the  one 
dimensional  Newton  algorithm  described  in  Section  5  of  this  report. 

7.3  Numerical  result.  We  now  consider  a  perforated  square  thin  plate  (width  «  1.), 
subjected  to  an  uniform  traction  of  .52  per  unit  area  on  two  of  its  opposite  faces.  This 
plate  is  supposed  to  be  made  of  a  Tresca  material  with  p  “  1.5  and  k  »  1//2  . 

For  symmetry  reasons,  only  one  fourth  of  the  plate  is  considered.  On  this  fourth, 

126  nodes  are  used  for  approximating  the  velocity  field.  One  block-relaxation  iteration 
is  done  per  Uzawa  step,  and  the  parameter  R  and  the  weight  r(x)  are  respectively  given 
by  R  -  1  and  r(x)  *  |h|p_2,  H  being  the  deformation  rate  tensor  corresponding  to  a 
computation  done  on  the  same  geometry  but  with  P^(G)  »  p(k/2)Pjc|P  (compressible  Horton 
material  in  plane  strains). 

After  50  iterations,  the  error  ID(»n)  -  B"»  is  decreased  by  a  factor  of  10~4,  and 
the  total  dissipated  energy  rate  la  equal  to  -2.737  for  the  whole  plate.  The 
corresponding  velocitltes  are  indicated  on  Fig  7,2.  It  must  be  noticed  that,  due  to  the 
little  number  of  boundary  conditions  imposed  on  v,  this  case  is  particularly  unstable  for 
most  numerical  methods. 

8  POSSIBLE  EXTENSIONS  OF  THU  METHOD. 

Many  extensions  can  be  considered  for  the  numerical  method  described  in  this 
report.  For  example, 

(i)  different  finite  elements  can  be  considered  in  the  approximation  Kj,  of  the  set 
of  kinematically  admissible  incompressible  velocity  fields.  Any  finite  element  which  is 
used  with  some  success  in  the  approximation  of  the  Stokes  problem  can  be  employed  here. 
Nevertheless,  if  the  gradients  of  the  elements  of  are  not  piecewise  constant,  a 

numerical  integration  rule  will  be  necessary  to  compute  the  dissipation  F(C),  which 
leads  to  an  additional  truncation  error  and  which  slightly  complicates  Problem  (3.6). 
Moreover,  the  space  Yh,  which  is  then  the  space  of  functions  which  are  characterized  by 
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their  values  at  the  integration  points,  must  be  sufficiently  large  to  contain,  within  an 
isomorphism,  the  image  of  by  the  operator  0< * ) » 

(ii)  an  inertia  term  can  be  added  in  the  formulation  of  the  virtual  work  theorem. 
Through  an  implicit  time  discretization,  the  resulting  problem  will  then  reduce  to  a 
sequence  of  augmented  lagrangian  problems  (3.6)  (one  per  time  step),  the  functional  G 
being  now  replaced  by 

G  (w)  -  ~l  f«w  dx  -  J  g*w  da  +  .■—  |  p  |*  -  v  |2  dx. 

a  r2  n  n 

Each  problem  (3.6)  can  still  be  solved  by  Algorithm  (3.11)  -  (3.14).  Problem  (3.12)  will 
again  correspond  to  a  linear  Stokes  type  problem,  associated  to  fixed,  symmetric,  positive 
definite  finite  element  matrices.  Problem  (3.13)  remains  unchanged; 

(iil)  a  convection  term  p(v*V)v  can  also  be  added  in  the  formulation  of  the  virtual 
work  theorem.  Since  the  operator  in  v  will  no  longer  be  self-adjoint,  no  augmented 
lagrangian  iR  can  then  be  introduced.  Nevertheless,  Algorithm  (3.11)  -  (3.14)  is  still 
applicable  there  ( FORTIN-GLOWINSKI  [1982,  p  71].  Problem  (3.13)  is  unchanged,  and  (3.12) 
becomes 

]  r(x)(R(B(v)-»)-X)-D(*)dx  +  I  p(v«V)vm  dx  -  j  f*m  dx  +  J  g>w  da,  V  <r  e  K  . 

n  n  n  r2 

For  small  convection  terms,  v  can  be  replaced  in  the  convection  term  by  the  solution  v" 
at  the  previous  iterate,  and  (3.12)  then  reduces  to  an  ordinary  Stokes  problem.  For  large 
convection  terms,  one  can  use  optimal  control  techniques  ( [GLOWINSKI-LE  TALLEC  [1983)). 

All  this  is  described  in  details  by  TANGUY  [1983]  which  uses  augmented  lagrangian 
techniques  in  a  very  similar  situation; 

(iv)  finally,  our  problem  can  be  coupled  to  an  heat  diffusion  problem  if  we  suppose, 
for  example,  that  the  internal  dissipation  potential  P^x,  G)  is  a  function  of  the 
temperature  T  at  x.  If  convection  phenomenon  are  not  dominant,  temperatures  and 
velocities  can  be  efficiently  computed  hy  block-relaxation:  assuming  the  velocity  to  be 
given,  one  computes  the  temperature  by  solving  the  energy  equation;  then,  assuming  the 


r 


temperature  to  be  given,  the  velocity  is  determined  by  solving  (3.6),  the  process  being 
repeated  until  convergence.  Observe  that,  despite  a  possible  change  of  the  temperature 
field  between  two  successive  resolutions  of  (3.6),  the  finite  element  matrices  do  not  have 
to  be  changed  because  the  temperature  is  only  a  parameter  in  the  local  problems  (4.1). 

This  results  in  considerable  economy  in  computer  running  time. 
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